Global Flight Route Density Mapping

🎯 What You’ll Create

By the end of this tutorial, you’ll have a beautiful, professional map showing where airplanes fly most frequently around the world. Think of it like a heat map - brighter colors show busier flight routes!

Perfect for: Beginners in data visualization, geography students, aviation enthusiasts, or anyone wanting to learn R programming.

📋 Overview

What This Tutorial Does (In Simple Terms)

1. Downloads real flight data from the internet (airports and routes)
2. Processes the data to connect airports with flight paths
3. Creates a hexagon grid over the world map (like a honeycomb)
4. Counts how many flights pass through each hexagon
5. Makes a colorful map showing flight density

Types of Routes We’ll Create

Straight Routes 🔹: Direct lines between airports (like drawing with a ruler on a flat map)

Curved Routes 🌐: Realistic paths following Earth’s spherical surface (how planes actually fly!)

NoteWhy Curved Routes Matter

Imagine you’re flying from New York to Tokyo. On a flat map, the shortest path looks like it goes straight across the Pacific Ocean. But Earth is round! The actual shortest path curves up through Alaska. That’s what geodetic (curved) routes show.

---

Install Required R Packages

Run this code once to install everything:

# Install all required packages
install.packages(c(
  "readr", "dplyr", "tidyr", "stringr", "sf", "ggplot2", 
  "viridis", "scales", "rnaturalearth", 
  "rnaturalearthdata", "geosphere", "tmap", "purrr"
))

TipInstallation Time

This may take 5-10 minutes depending on your internet speed. Be patient!

Load Libraries

After installing packages once, you need to “load” them each time you start R:

# Load all libraries we'll use
library(readr)
library(dplyr)
library(tidyr)
library(stringr)    
library(sf)      
library(ggplot2)       
library(viridis)      
library(scales)      
library(rnaturalearth) 
library(geosphere)     
library(tmap)
library(purrr)     

cat("✓ All libraries loaded successfully!\n")

✓ All libraries loaded successfully!

---

📂 Step 1: Download and Prepare Airport Data

Create Folder Structure

# Create two folders: one for outputs, one for raw data
folders <- c("output", "data")
sapply(folders, dir.create, showWarnings = FALSE, recursive = TRUE)

output   data 
  TRUE   TRUE 

cat("✓ Folders created successfully!\n")

✓ Folders created successfully!

Download Airport Data

Download the airport location data from the jpatokal/openflights repository:

# Check if we already have the airports file
if (!file.exists("data/airports.dat")) {
  cat("📥 Downloading airports.dat...\n")
  
  download.file(
    url = "https://raw.githubusercontent.com/jpatokal/openflights/master/data/airports.dat",
    destfile = "data/airports.dat",
    mode = "wb"
  )
  cat("✓ Airports downloaded!\n")
} else {
  cat("✓ Airports file already exists. Skipping download.\n")
}

📥 Downloading airports.dat...
✓ Airports downloaded!

Download Routes Data

Download the routes data from the jpatokal/openflights repository:

# Check if we already have the routes file
if (!file.exists("data/routes.dat")) {
  cat("📥 Downloading routes.dat...\n")
  
  download.file(
    url = "https://raw.githubusercontent.com/jpatokal/openflights/master/data/routes.dat",
    destfile = "data/routes.dat",
    mode = "wb"
  )
  cat("✓ Routes downloaded!\n")
} else {
  cat("✓ Routes file already exists. Skipping download.\n")
}

📥 Downloading routes.dat...
✓ Routes downloaded!

Load World Boundaries

Load world boundaries from the rnaturalearth library:

# Get world boundaries
world <- ne_countries(scale = "medium", returnclass = "sf") %>%
  st_transform(8857) %>%  # Transform to Equal Earth projection
  st_make_valid()         # Ensure geometries are valid

cat("✓ Loaded world boundaries\n")

✓ Loaded world boundaries

NoteCoordinate Reference System (CRS)

We use EPSG:8857 (Equal Earth projection) which is excellent for global maps. It balances shape, area, and distance distortion across the entire globe.

Load Routes Data

# Define column names for routes
routes_column <- c(
  "Airline", "Airline_ID", "Source_airport", "Source_airport_ID",
  "Destination_airport", "Destination_airport_ID", 
  "Codeshare", "Stops", "Equipment"
)

# Read the routes file
routes <- read_delim(
  "data/routes.dat", 
  delim = ",",
  col_names = routes_column,
  show_col_types = FALSE
)

# Save a copy as CSV
write_csv(routes, "output/routes.csv")

cat("✓ Loaded", nrow(routes), "routes from the database\n")

✓ Loaded 67663 routes from the database

Load Airports Data

# Define column names for airports
airports_column <- c(
  "Airport_ID", "Name", "City", "Country", "IATA", "ICAO",
  "Latitude", "Longitude", "Altitude", "Timezone", "DST",
  "Tz_database_timezone", "Type", "Source"
)

# Read the airports file
airports <- read_delim(
  "data/airports.dat", 
  delim = ",",
  col_names = airports_column,
  show_col_types = FALSE
)

# Save a copy as CSV
write_csv(airports, "output/airports.csv")

cat("✓ Loaded", nrow(airports), "airports from the database\n")

✓ Loaded 7698 airports from the database

cat("\n📍 Sample airports:\n")

📍 Sample airports:

print(head(airports %>% select(Name, City, Country, IATA, Latitude, Longitude), 5))

# A tibble: 5 × 6
  Name                                    City  Country IATA  Latitude Longitude
  <chr>                                   <chr> <chr>   <chr>    <dbl>     <dbl>
1 Goroka Airport                          Goro… Papua … GKA      -6.08      145.
2 Madang Airport                          Mada… Papua … MAG      -5.21      146.
3 Mount Hagen Kagamuga Airport            Moun… Papua … HGU      -5.83      144.
4 Nadzab Airport                          Nadz… Papua … LAE      -6.57      147.
5 Port Moresby Jacksons International Ai… Port… Papua … POM      -9.44      147.

Convert Airports to Spatial Format

Convert the airports dataframe to a spatial object:

ImportantData Requirements

The Latitude and Longitude fields must not contain missing values for proper conversion to spatial format.

# Convert to spatial object
airports_sf <- airports %>%
  mutate(
    Latitude = as.numeric(Latitude),
    Longitude = as.numeric(Longitude)
  ) %>%
  filter(!is.na(Latitude), !is.na(Longitude)) %>%
  st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326, remove = FALSE)

# Transform to Equal Earth projection
airports_sf <- st_transform(airports_sf, crs = 8857)

# Save as shapefile
st_write(airports_sf, "output/airports.shp", delete_dsn = TRUE, quiet = TRUE)

cat("✓ Created spatial object with", nrow(airports_sf), "airports\n")

✓ Created spatial object with 7698 airports

Visualize Airport Locations

# Create airport locations map
airports_map <- tm_shape(world) + 
  tm_polygons(
    col = "gray20", 
    border.col = NA,
    lwd = 0
  ) +  
  tm_shape(airports_sf) + 
  tm_dots(
    col = "#00B2D1",
    size = 0.1,
    alpha = 0.6,
    border.lwd = 0
  ) + 
  tm_layout(
    title = "GLOBAL AIRPORT LOCATIONS",
    title.size = 1.2,
    title.color = "gray15",
    title.position = c("center", "top"),
    title.fontface = "bold",
    bg.color = "gray90",
    frame = FALSE,
    legend.position = c("left", "bottom"),
    legend.bg.color = "#16213e",
    legend.bg.alpha = 0.85,
    legend.frame.color = "#0f3460",
    legend.frame.lwd = 2,
    legend.title.size = 1.1,
    legend.title.color = "#ffffff",
    legend.title.fontface = "bold",
    legend.text.size = 0.85,
    legend.text.color = "#e8e8e8",
    inner.margins = c(0.02, 0.02, 0.08, 0.02),
    outer.margins = c(0, 0, 0, 0),
    asp = 1.8,
    fontfamily = "sans"
  ) +
  tm_credits(
    text = "Source: jpatokal/openflights | Year: 2014 | Total Airports: 7,698", 
    position = c("center", "bottom"), 
    col = "gray40",
    size = 0.6
  )

# Save the map
tmap_save(airports_map, "output/01.airports_map.png", width = 18, height = 10,  dpi = 400)

# Display the map
airports_map

---

🛫 Step 2: Create Route Lines

Clean Routes Data

Clean the routes data before joining with airports:

# Clean the routes data
routes_clean <- routes %>%
  mutate(
    Source_airport = str_trim(Source_airport),
    Destination_airport = str_trim(Destination_airport)
  ) %>%
  filter(
    !is.na(Source_airport), !is.na(Destination_airport),
    Source_airport != "", Destination_airport != "",
    str_detect(Source_airport, "^[A-Z]{3}$"),
    str_detect(Destination_airport, "^[A-Z]{3}$")
  )

cat("✓", nrow(routes_clean), "valid routes (from", nrow(routes), "total)\n")

✓ 67663 valid routes (from 67663 total)

Join Airport Coordinates to Routes

# Create lookup table for airport coordinates
airports_lookup <- airports %>%
  mutate(
    IATA = str_trim(IATA),
    Latitude = as.numeric(Latitude),
    Longitude = as.numeric(Longitude)
  ) %>%
  filter(
    !is.na(IATA), IATA != "",
    str_detect(IATA, "^[A-Z]{3}$"),
    !is.na(Latitude), !is.na(Longitude)
  ) %>%
  distinct(IATA, .keep_all = TRUE) %>%
  select(IATA, Latitude, Longitude)

# Join coordinates to routes
routes_joined <- routes_clean %>%
  left_join(
    airports_lookup %>% rename(src_lat = Latitude, src_lon = Longitude),
    by = c("Source_airport" = "IATA")
  ) %>%
  left_join(
    airports_lookup %>% rename(dst_lat = Latitude, dst_lon = Longitude),
    by = c("Destination_airport" = "IATA")
  ) %>%
  filter(
    !is.na(src_lon), !is.na(src_lat),
    !is.na(dst_lon), !is.na(dst_lat)
  )

write_csv(routes_joined, "output/routes_inc_locations.csv")
cat("✓ Saved routes with coordinates:", nrow(routes_joined), "routes\n")

✓ Saved routes with coordinates: 66934 routes

Count Airport Traffic

Count how many routes use each airport as source or destination:

# Count source and destination airports
src_count <- routes_joined %>% count(Source_airport) %>% rename(IATA = Source_airport, src_count = n)

dst_count <- routes_joined %>% count(Destination_airport) %>% rename(IATA = Destination_airport, dst_count = n)

# Combine counts
airports_count <- full_join(src_count, dst_count, by = "IATA") %>%
  mutate(
    src_count = replace_na(src_count, 0),
    dst_count = replace_na(dst_count, 0),
    total = src_count + dst_count
  )

# Join with coordinates
airports_count <- left_join(airports_lookup, airports_count, by = "IATA")

# Convert to spatial object
airports_count_sf <- airports_count %>%
  filter(!is.na(Latitude), !is.na(Longitude)) %>%
  st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326, remove = FALSE) %>%
  st_transform(crs = 8857)

cat("✓ Calculated airport traffic counts\n")

✓ Calculated airport traffic counts

cat("📊 Top 5 busiest airports by route count:\n")

📊 Top 5 busiest airports by route count:

print(airports_count %>% arrange(desc(total)) %>% head(5) %>% select(IATA, total))

# A tibble: 5 × 2
  IATA  total
  <chr> <int>
1 ATL    1826
2 ORD    1108
3 LHR    1051
4 PEK    1050
5 CDG    1041

Visualize Airport Traffic

# Prepare the map
world_busiest_airports_map <- tm_shape(world, bbox = st_bbox(c(xmin = -180, xmax = 180, ymin = -90, ymax = 90), crs = st_crs(4326))) + 
  tm_polygons(
    col = "#525252", 
    border.col = "#252525",
    lwd = 0.1
  ) + 
  tm_shape(airports_count_sf) + 
  tm_dots(
    col = "total",        
    size = "total",         
    palette = "viridis",     
    style = "quantile",
    size.scale = tm_scale_continuous(values.scale = 1000),
    border.lwd = 0,                 
    alpha = 0.6,
    legend.reverse = TRUE,
    legend.col.show = TRUE,
    legend.size.show = FALSE,
    col.legend = tm_legend(title = "Total Flights")
  ) + 
  tm_layout(
    title = "WORLD'S BUSIEST AIRPORTS",
    title.size = 1.2,
    title.color = "#ffffff",
    title.position = c("center", "top"),
    title.fontface = "bold",
    bg.color = "#252525",
    frame = FALSE,
    legend.position = c("left", "bottom"),
    legend.bg.color = "#16213e",
    legend.bg.alpha = 0.85,
    legend.frame.color = "#0f3460",
    legend.frame.lwd = 2,
    legend.title.size = 1.1,
    legend.title.color = "#ffffff",
    legend.title.fontface = "bold",
    legend.text.size = 0.85,
    legend.text.color = "#e8e8e8",
    inner.margins = c(0.02, 0.02, 0.08, 0.02),
    outer.margins = c(0, 0, 0, 0),
    asp = 1.8,
    fontfamily = "sans"
  ) +
  tm_credits(
    text = "Source: jpatokal/openflights | Year: 2014", 
    position = c("center", "bottom"), 
    col = "gray40",
    size = 0.6
  )

# Save the map
tmap_save(world_busiest_airports_map, "output/02.world_busiest_airports_map.png", width = 18, height = 10,  dpi = 400)

# Display the map
world_busiest_airports_map

Create Straight (Euclidean) Route Lines

Create direct line connections between airports:

# Create straight line geometries
routes_lines_sf <- routes_joined %>%
  mutate(
    geometry = mapply(
      function(x1, y1, x2, y2) {
        st_linestring(matrix(c(x1, y1, x2, y2), ncol = 2, byrow = TRUE))
      },
      src_lon, src_lat, dst_lon, dst_lat,
      SIMPLIFY = FALSE
    ) %>% st_sfc(crs = 4326)
  ) %>%
  st_as_sf() %>%
  st_transform(crs = 8857)

# Save as shapefile
st_write(routes_lines_sf, "output/routes_lines.shp", delete_dsn = TRUE, quiet = TRUE)

cat("✓ Created", nrow(routes_lines_sf), "straight route lines\n")

✓ Created 66934 straight route lines

Visualize Straight Routes

# Create map with tmap
routes_lines_map <- tm_shape(world) + 
  tm_polygons(
    col = "#525252", 
    border.col = "#252525",
    lwd = 0.1
  ) +  
  tm_shape(routes_lines_sf) + 
  tm_lines(
    col = "#ffffff",
    lwd = 0.1, 
    alpha = 0.08
  ) +
  tm_shape(airports_count_sf) + 
  tm_dots(
    col = "total",        
    size = "total",         
    palette = "viridis",     
    style = "quantile",
    size.scale = tm_scale_continuous(values.scale = 10),
    border.lwd = 0,                 
    alpha = 0.6,
    legend.reverse = TRUE,
    legend.col.show = TRUE,
    legend.size.show = FALSE,
    col.legend = tm_legend(title = "Total Flights")
  ) +   
  tm_layout(
    title = "GLOBAL AIRPLANE ROUTES (EUCLIDEAN DISTANCE)",
    title.size = 1.2,
    title.color = "#ffffff",
    title.position = c("center", "top"),
    title.fontface = "bold",
    bg.color = "#252525",
    frame = FALSE,
    legend.position = c("left", "bottom"),
    legend.bg.color = "#16213e",
    legend.bg.alpha = 0.85,
    legend.frame.color = "#0f3460",
    legend.frame.lwd = 2,
    legend.title.size = 1.1,
    legend.title.color = "#ffffff",
    legend.title.fontface = "bold",
    legend.text.size = 0.85,
    legend.text.color = "#e8e8e8",
    inner.margins = c(0.02, 0.02, 0.08, 0.02),
    outer.margins = c(0, 0, 0, 0),
    asp = 1.8,
    fontfamily = "sans"
  ) +
  tm_credits(
    text = "Source: jpatokal/openflights | Year: 2014", 
    position = c("center", "bottom"), 
    col = "gray80",
    size = 0.6
  )

# Save the map
tmap_save(routes_lines_map, "output/03.routes_euclidean_map.png", width = 18, height = 10,  dpi = 400)

# Display the map
routes_lines_map

WarningWhy Straight Lines Are Misleading

While straight lines look simple on a flat map, they don’t represent how planes actually fly. Earth’s curvature means the shortest distance between two points is a curve, not a straight line!

Create Curved (Geodetic) Route Lines

Create realistic curved routes following Earth’s surface:

# Convert linestring to curves
routes_geodetic_sf <- routes_joined %>%
  mutate(
    geometry = pmap(
      list(src_lon, src_lat, dst_lon, dst_lat),
      function(lon1, lat1, lon2, lat2) {
        # Calculate great circle intermediate points
        gc_points <- gcIntermediate(
          c(lon1, lat1), 
          c(lon2, lat2), 
          n = 100,  # Number of intermediate points
          addStartEnd = TRUE
        )
        st_linestring(gc_points)
      }
    ) %>% st_sfc(crs = 4326)
  ) %>%
  st_as_sf() %>%
  st_transform(crs = 8857)

# Save as shapefile
st_write(routes_geodetic_sf, "output/routes_geodetic_lines.shp", delete_dsn = TRUE, quiet = TRUE)

cat("✓ Created", nrow(routes_geodetic_sf), "geodetic routes\n")

✓ Created 66934 geodetic routes

TipHow Geodetic Lines Work

The gcIntermediate() function calculates the shortest path between two points on a sphere (Earth). By adding 100 intermediate points, we get a smooth curve that accurately represents the flight path.

Visualize Geodetic Routes

# Create map with tmap
routes_geodetic_map <- tm_shape(world) + 
  tm_polygons(
    col = "#525252", 
    border.col = "#252525", 
    lwd = 0
  ) +  
  tm_shape(routes_geodetic_sf) + 
  tm_lines(
    col = "#ffffff", 
    lwd = 0.1, 
    alpha = 0.08
  ) + 
  tm_shape(airports_count_sf) + 
  tm_dots(
    col = "total",        
    size = "total",         
    palette = "viridis",     
    style = "quantile",
    size.scale = tm_scale_continuous(values.scale = 10),
    border.lwd = 0,                 
    alpha = 0.6,
    legend.reverse = TRUE,
    legend.col.show = TRUE,
    legend.size.show = FALSE,
    col.legend = tm_legend(title = "Total Flights")
  ) + 
  tm_layout(
    title = "GLOBAL AIRPLANE ROUTES (GEODESIC DISTANCE)",
    title.size = 0.9,
    title.color = "#ffffff",
    title.position = c("center", "top"),
    title.fontface = "bold",
    bg.color = "#252525",
    frame = FALSE,
    legend.position = c("left", "bottom"),
    legend.bg.color = "#16213e",
    legend.bg.alpha = 0.85,
    legend.frame.color = "#0f3460",
    legend.frame.lwd = 2,
    legend.title.size = 1.1,
    legend.title.color = "#ffffff",
    legend.title.fontface = "bold",
    legend.text.size = 0.85,
    legend.text.color = "#e8e8e8",
    inner.margins = c(0.02, 0.02, 0.08, 0.02),
    outer.margins = c(0, 0, 0, 0),
    asp = 1.8,
    fontfamily = "sans"
  ) +
  tm_credits(
    text = "Source: jpatokal/openflights | Year: 2014", 
    position = c("center", "bottom"), 
    col = "#a0a3ad", 
    size = 0.7
  )

# Save the map
tmap_save(routes_geodetic_map, "output/04.routes_geodetic_map.png", width = 18, height = 10,  dpi = 400)

# Display the map
routes_geodetic_map

---

🔶 Step 3: Create Hexagonal Grid

Generate Hexagon Grid

• In this study, a hexagonal grid overlay was created for the world.
• The grid has a cell size of 300 km, specified as cellsize = 3e5 (300,000 meters or 300 km).
• This grid provides a spatial framework for analyzing the global distribution of flight routes and other spatial data.
• This type of hexagonal or square grid overlay is particularly useful when the data density is high, and traditional visualization techniques are unable to accurately depict the true variation in the data. By dividing the space into uniform grid cells, it becomes easier to identify patterns and trends, especially in areas with complex or dense data.

# Create hexagonal grid
hexgrid <- st_make_grid(
  world,
  cellsize = 3e5,
  what = "polygons",
  square = FALSE   # FALSE creates hexagons
) %>%
  st_as_sf() %>%
  mutate(id = row_number())

cat("✓ Created hexagon grid with", nrow(hexgrid), "hexagons\n")

✓ Created hexagon grid with 7557 hexagons

NoteWhy Hexagons?

Hexagons are ideal for spatial analysis because:

• Each hexagon touches 6 neighbors (more than squares)
• They minimize edge effects and sampling bias
• They’re visually appealing and easy to interpret

Filter to Land Hexagons

Keep only hexagons that intersect with land:

# Union all countries into single geometry
land_union <- st_union(world)

# Keep hexagons that intersect land
hexgrid_world <- hexgrid[lengths(st_intersects(hexgrid, land_union)) > 0, ]

cat("✓ Filtered to", nrow(hexgrid_world), "hexagons over land\n")

✓ Filtered to 2782 hexagons over land

# Save hexgrid
st_write(hexgrid_world, "output/world_hexgrid.shp", delete_layer = TRUE, quiet = TRUE)

---

🔢 Step 4: Count Routes per Hexagon

Calculate Route-Hexagon Intersections

Count how many routes pass through each hexagon:

cat("🔍 Calculating route-hexagon intersections...\n")

🔍 Calculating route-hexagon intersections...

cat("⏱️ This may take 1-2 minutes for", nrow(routes_geodetic_sf), "routes\n\n")

⏱️ This may take 1-2 minutes for 66934 routes

# Find which routes intersect each hexagon
intersections_geodetic <- st_intersects(hexgrid_world, routes_geodetic_sf)

cat("✓ Intersection calculation complete\n")

✓ Intersection calculation complete

Add Route Counts to Hexagons

# Count routes per hexagon
hexgrid_world_geodetic_cnt <- hexgrid_world %>% mutate(nm_rts = lengths(intersections_geodetic))

# Save result
st_write(hexgrid_world_geodetic_cnt, "output/world_hexgrid_geodetic_count.shp", delete_layer = TRUE, quiet = TRUE)

cat("✓ Added route counts to hexagons\n\n")

✓ Added route counts to hexagons

# Summary statistics
cat("📊 Route Count Statistics:\n")

📊 Route Count Statistics:

summary_stats <- summary(hexgrid_world_geodetic_cnt$nm_rts)
print(summary_stats)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    0.0    40.0    82.0   207.1   205.8  4954.0 

# Additional insights
cat("\n📍 Key Insights:\n")

📍 Key Insights:

cat("• Hexagons with routes:", sum(hexgrid_world_geodetic_cnt$nm_rts > 0), "\n")

• Hexagons with routes: 2497 

cat("• Hexagons with 100+ routes:", sum(hexgrid_world_geodetic_cnt$nm_rts >= 100), "\n")

• Hexagons with 100+ routes: 1232 

cat("• Maximum routes in a hexagon:", max(hexgrid_world_geodetic_cnt$nm_rts), "\n")

• Maximum routes in a hexagon: 4954 

---

🎨 Step 5: Create Final Visualization

Create Flight Density Heat Map

# Create the map using tmap
map_geodetic <- tm_shape(world) +
  tm_polygons(
    col = "#1a1a2e",
    border.col = "#0f3460",
    lwd = 0.3,
    alpha = 0.6
  ) +
  tm_shape(hexgrid_world_geodetic_cnt) +
  tm_polygons(
    col = "nm_rts",
    palette = c("#440154", "#31688e", "#35b779", "#fde724"), 
    style = "fixed",
    alpha = 0.9,
    border.col = "#0a0e27",
    lwd = 2,
    title = "Flight Routes",
    breaks = c(0, 10, 50, 100, 250, 500, 1000, 2500, max(hexgrid_world_geodetic_cnt$nm_rts)),
    labels = c("0-10", "10-50", "50-100", "100-250", "250-500", "500-1000", "1000-2500", paste0(">2500")),
    legend.reverse = TRUE,
    legend.format = list(text.separator = " ")
  ) +
  tm_layout(
    title = "Global Flight Route Density",
    title.size = 1.4,
    title.color = "#ffffff",
    title.position = c("center", "top"),
    title.fontface = "bold",
    bg.color = "#0a0e27",
    frame = FALSE,
    legend.position = c("left", "bottom"),
    legend.bg.color = "#16213e",
    legend.bg.alpha = 0.85,
    legend.frame.color = "#0f3460",
    legend.frame.lwd = 2,
    legend.title.size = 1.1,
    legend.title.color = "#ffffff",
    legend.title.fontface = "bold",
    legend.text.size = 0.85,
    legend.text.color = "#e8e8e8",
    inner.margins = c(0.02, 0.02, 0.08, 0.02),
    outer.margins = c(0, 0, 0, 0),
    asp = 1.8,
    fontfamily = "sans"
  ) +
  tm_credits(
    text = "Data: OpenFlights (2014) | Year: 2014", 
    position = c("center", "bottom"), 
    col = "#94a3b8", 
    size = 0.75,
    fontface = "italic"
  ) +
  tm_compass(
    type = "arrow",
    position = c("right", "top"),
    size = 2,
    color.dark = "#fde724",
    color.light = "#ffffff",
    text.size = 0.6,
    text.color = "#ffffff"
  ) +
  tm_scale_bar(
    position = c("right", "bottom"),
    text.size = 0.6,
    text.color = "#ffffff",
    color.dark = "#ffffff",
    color.light = "#94a3b8",
    breaks = c(0, 2000, 4000, 6000)
  )

# Save the enhanced map
tmap_save(map_geodetic,"output/05.Hexagon_Grid_flight_density_map.png",width = 18,height = 10,dpi = 400)

# Display the map
map_geodetic

TipColor Scale Interpretation

• Dark purple/blue (0-50 routes): Remote or less-connected regions
• Pink/orange (50-250 routes): Moderate air traffic
• Bright yellow (250+ routes): Major aviation corridors and hubs

---

🎉 Success!

Congratulations! You’ve successfully created a professional flight route density visualization!

What You’ve Learned

✅ How to download and process geospatial data
✅ The difference between Euclidean and geodetic distances
✅ How to create and work with hexagonal grids
✅ How to perform spatial intersections
✅ How to create professional data visualizations with tmap

Key Insights from Your Map

Bright colors (yellow/pink) indicate:

• High concentration of flight routes
• Major aviation hubs (e.g., Europe, North America, East Asia)
• Well-connected regions with frequent air traffic

Dark colors (purple/blue) indicate:

• Few flight routes
• Remote or less connected regions
• Lower air traffic density

Next Steps

Want to take this further? Try:

1. Filter by airline or region: Analyze specific airline networks
2. Time-based analysis: Compare route patterns across years
3. Network analysis: Calculate connectivity metrics between airports
4. Interactive maps: Use leaflet to create web-based interactive versions

---

📚 Additional Resources

• sf package documentation
• tmap documentation
• OpenFlights database
• Spatial data science book

---

📝 Session Information

sessionInfo()

R version 4.5.2 (2025-10-31 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 11 x64 (build 26200)

Matrix products: default
  LAPACK version 3.12.1

locale:
[1] LC_COLLATE=English_United States.utf8 
[2] LC_CTYPE=English_United States.utf8   
[3] LC_MONETARY=English_United States.utf8
[4] LC_NUMERIC=C                          
[5] LC_TIME=English_United States.utf8    

time zone: Europe/Berlin
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] purrr_1.2.0         tmap_4.2            geosphere_1.5-20   
 [4] rnaturalearth_1.1.0 scales_1.4.0        viridis_0.6.5      
 [7] viridisLite_0.4.2   ggplot2_4.0.1       sf_1.0-23          
[10] stringr_1.6.0       tidyr_1.3.1         dplyr_1.1.4        
[13] readr_2.1.6        

loaded via a namespace (and not attached):
 [1] tidyselect_1.2.1        farver_2.1.2            S7_0.2.1               
 [4] fastmap_1.2.0           leaflegend_1.2.1        leaflet_2.2.3          
 [7] XML_3.99-0.20           digest_0.6.37           lifecycle_1.0.4        
[10] terra_1.8-86            magrittr_2.0.4          compiler_4.5.2         
[13] rlang_1.1.6             tools_4.5.2             utf8_1.2.6             
[16] yaml_2.3.11             data.table_1.17.8       knitr_1.50             
[19] htmlwidgets_1.6.4       bit_4.6.0               sp_2.2-0               
[22] classInt_0.4-11         RColorBrewer_1.1-3      abind_1.4-8            
[25] KernSmooth_2.23-26      withr_3.0.2             leafsync_0.1.0         
[28] grid_4.5.2              cols4all_0.10           e1071_1.7-16           
[31] leafem_0.2.5            colorspace_2.1-2        spacesXYZ_1.6-0        
[34] cli_3.6.5               crayon_1.5.3            rmarkdown_2.30         
[37] generics_0.1.4          tzdb_0.5.0              tmaptools_3.3          
[40] DBI_1.2.3               proxy_0.4-27            stars_0.7-0            
[43] parallel_4.5.2          s2_1.1.9                base64enc_0.1-3        
[46] vctrs_0.6.5             jsonlite_2.0.0          hms_1.1.4              
[49] bit64_4.6.0-1           crosstalk_1.2.2         units_1.0-0            
[52] maptiles_0.10.0         glue_1.8.0              lwgeom_0.2-14          
[55] codetools_0.2-20        stringi_1.8.7           gtable_0.3.6           
[58] raster_3.6-32           tibble_3.3.0            logger_0.4.1           
[61] pillar_1.11.1           htmltools_0.5.8.1       rnaturalearthdata_1.0.0
[64] R6_2.6.1                wk_0.9.4                microbenchmark_1.5.0   
[67] vroom_1.6.6             evaluate_1.0.5          lattice_0.22-7         
[70] png_0.1-8               class_7.3-23            Rcpp_1.1.0             
[73] gridExtra_2.3           xfun_0.54               pkgconfig_2.0.3        

---

💾 Files Created

This tutorial creates the following files in the output folder:

File	Description
routes.csv	Raw routes data
airports.csv	Raw airports data
routes_inc_locations.csv	Routes with coordinate information
airports.shp	Airports as spatial points
routes_lines.shp	Straight (Euclidean) route lines
routes_geodetic_lines.shp	Curved (geodetic) route lines
world_hexgrid.shp	Hexagonal grid over land
world_hexgrid_geodetic_count.shp	Hexagons with route counts
01.airports_map.png	Map showing all airport locations
02.world_busiest_airports_map.png	Map of busiest airports by route count
03.routes_euclidean_map.png	Map with straight routes and airports
04.routes_geodetic_map.png	Map with curved routes and airports
05.Hexagon_Grid_flight_density_map.png	Final heat map visualization

---

End of Tutorial ✈️🌍